function dz=foss(t,z)
 
global parms
% Note that parms is a 1x14 vector with elements:
%
% parms = [delta A Sbar alpha0 alpha1 alpha2 alpha3 Gamma1 alpha psi beta gamma Q0 popgr];
%           1    2   3      4      5    6     7       8      9    10  11    12  13  14  
%
% 
delta = parms(1);
A = parms(2);
Sbar = parms(3);
alpha0 = parms(4);
alpha1 = parms(5);
alpha2 = parms(6);
alpha3 = parms(7);
beta = 1/parms(11)-1;
gamma = 12;
popgr = parms(14);
 
Q = exp(popgr*t);
 
dz=zeros(5,1);
 
k=z(1);
S=z(2);
N=z(3);
lambda=z(4);
sigma=z(5);
nu=z(6);
 
c = lambda^(-1/gamma);
g = alpha0 + alpha1/(Sbar-alpha2/(alpha3+N)-S);
gdS = alpha1*(alpha3+N)^2/((Sbar-S)*(alpha3+N)-alpha2)^2;
gdN = -alpha1*alpha2/((Sbar-S)*(alpha3+N)-alpha2)^2;
gd2N = 2*alpha1*alpha2*(Sbar-S)/((Sbar-S)*(alpha3+N)-alpha2)^3;
gdSN = -2*alpha1*alpha2*(alpha3+N)/((Sbar-S)*(alpha3+N)-alpha2)^3;
 
%
% Now solve separately for cases where n = 0 or n > 0
%

if lambda > nu
	n = 0;
    dz(6) = beta*nu+gdN*A*lambda.*k;
      
else

	n=(gdN*k*(delta+g*A-A+sigma*A*Q/lambda)+c*gdN-gdSN*Q*A*k^2-...
        		popgr*sigma.*Q/lambda)./(gd2N*k-2*gdN);
            dz(6) = lambda*(beta+delta-A+g*A)-sigma*Q*A;
end

i=A*k.*(1-g)-c-n;
 
dz(1) = i-delta*k;
dz(2) = A*Q*k;
dz(3) = n;
dz(4) = lambda*(beta+delta-A+g*A)-sigma*Q*A;
dz(5) = beta*sigma+gdS*A*lambda.*k;
